Identification of metabolizing enzyme genes associated with xenobiotics and odorants in the predatory stink bug Arma custos based on transcriptome analysis

The predatory stink bug, Arma custos, is a highly effective beneficial predator of crop pests. The lack of gene information related to xenobiotic detoxification and odorant degrading enzymes in the predator stink bugs to date has limited our ability for more in-depth studies of biological control. Hence, we conducted de novo assembly of the A. custos transcriptome from guts, antennae, and other tiussue samples of 5th instar larvae using Illumina sequencing technology. A total of 91, 50 and 23 genes of cytochrome P450 monooxygenases (CYPs), carboxyl/choline esterases (CCEs) and glutathione S-transferases (GSTs) genes were identified, respectively. Gene expansions of CYP3 and CYP4 clans and the hormone and pheromone processing CCE class were found in A. custos. Analysis of tissue-specific expression patterns showed that 37 CYPs, 14 CCEs and 8 GSTs were enriched in guts, and 6 CYPs, 5 CCEs and 2 GSTs were up-regulated in antennae, suggesting their potential roles on xenobiotics detoxification and ordorant degradation. Gene information data presented here could be useful for a deeper understanding of the ecology, physiology and behavior of this beneficial species and could be helpful to improve their bio-control efficiency.


Introduction
Arma custos Fallou (Hemiptera: Pentatomidae) is a predaceous insect species that is found throughout China, Mongolia, the Korean peninsula, and other Eastern Asian regions [1,2]. It can effectively control a wide range of agricultural and forest insect pests, such as those in the orders Lepidoptera, Coleoptera, Hymenoptera and Hemiptera [2]. For this reason, A. custos has been studied as a biological control agent in integrated pest management (IPM) program. Currently in China, A. custos can be sucessfully massly bred by a number of nature enemy factories. For example, in Guizhou province, more than six millions of A. custos were produced annually, and releasing A. custos has been proven to be especially effective in the biological control of tobacco pests, such as Spodoptera litura and Helicoverpa assultam, and of S. frugiperda in corn and sorghum fields [3].
The importance of this insect for biological control has stimulated several studies into its biology and ecology, including its taxonomy, morphology, reproduction, geographical notes, olfactory system, predator-prey interaction, feeding strategies, mitochondrial genome and microRNAs [2,[4][5][6][7][8][9][10][11][12]. However, the general knowledge about A. custos has important gaps. Currently releasing biological control program of this predatory stink bug is challenged by two issues. Firstly, the predatory bugs that are long-stem bred in factories often exhibit unsatisfactory performances in the establishment and control of the target pests in the fields. The olfactory system is critical for their poor adaptation to unfamiliar foods and environments. Wu et al. [6] identified chemosensory genes of A. custos, including odorant-binding proteins (OBPs), chemosensory proteins (CSPs) and chemosensory receptors. Odorant degrading enzymes (ODEs) could also be of major importance in olfaction dynamics, however, ODEs haven't been identified in A. custos. ODEs play pivotal roles in the inactive metabolism of exogenous odorants and in the recovery of sensitivity in the olfactory system to detect new odorants [13]. Until now, only one study was published with the genetic information about ODEs in a predaroty bug, Picromerus lewisi [14]. Further understanding the molecular mechanisms of the olfaction system including ODEs in A. custos could help us better answer these questions about field establishment of the predatory bugs.
What's more, biological control is usually accompanied by chemical control for practical applications in the field. These predatory stink bugs released in the wild fields are very difficult to avoid the exposure to insecticides. Therefore, the risk of insecticides on A. custos should not be ignored in IPM. Information of insecticide toxicities against the predators is helpful for exploiting insecticides which do not harm predators and subsequently improve the efficacy of IPM strategies. Our previous studies showed that imidacloprid, thiamethoxam and pyrethoids exhibited high toxicities against A. custos while sulfoxaflor, chlorantraniliprole and a couple of bioinsecticides had relatively low toxicitites against A. custos [3]. The genetic information about insecticide detoxification was revealed in a couple of predatory bugs, including Cyrtorhinus lividipennis, P. lewisi, and Eocanthecona furcellata [14][15][16][17]. However, genetic information on A. custos in response to chemical insecticides is lacking. Only two recent studies had assessed differential expression profiles of the detoxification-related genes in A. custos under exposure to Bacillus thuringiensis (Bt) toxins and Transgenic Bt cotton [18,19]. In-depth investigation of mechanisms of insecticide detoxification of the predators could help us to develop successful IPM strategies based on a better balance between chemical and biological control methods.
Insects rely mainly on three superfamilies of metabolizing enzymes to degrade environmental chemicals including odorants, insecticides and other xenobiotics: cytochrome P450 monooxygenases (CYPs), carboxyl/choline esterases (CCEs) and glutathione Stransferases (GSTs) [20]. However, the genetic information related to these metabolizing enzymes of A. custos is lacking. The aim of the work reported here was to uncover the genes of these three metabolizing enzyme families in A. custos by analyzing the transcriptomics and expression profile data. These results provide a comprehensive foundation for further studies of molecular mechanisms of processing environmental chemicals in the predatory stink bug, which could be helpful to develop better strategies for the use of these biological control agents.

Insects
The A. custos individuals used in the present study belonged to a colony reared in the natural enemy breeding center at Fenggang County of The Guizhou Provincial Tobacco Company Zunyi Branch, Zunyi, China. The bugs were reared continuously on Mythimna separata larvae that were fed with artificial diet mainly consisting of corn leaf powder. They were reared under 27 ± 1 • C, RH of 75 ± 5%, and a 16 h:8 h (L:D) photoperiod.

Sample preparation and RNA sequencing
Guts (AcG), antennae (AcA), salivary glands (AcSG), legs (AcL) and heads without antennae and salivary glands (AcH) from 100 nymphs at the fifth instar stage of A. custos were isolated for total RNA extraction using TRIzol reagent (Invitrogen, USA). Nymphs of A. custos were starved for at least 12 h prior to extraction of RNAs. Three biological replicates were perfomed on each tissue sample. In this study, RNA sequencing (RNA-seq) was carried out at Novogene Bioinformatics Technology Co., China. Using illumina NovaSeq 6000 (illumina, USA), the end reading of 150 bp pairing is generated. Each sample generated approximately 20 million sequence reads.
Raw data were firstly processed through in-house perl scripts by removing reads containing adapter or poly-N and low-quality reads. At the same time, Q20, Q30 and GC content of the clean data were calculated. The clean reads were then assembled into unigenes using the Trinity software (v2.6.6) [21].
Databases of the non-redundant protein sequences of National Center for Biotechnology Information (NCBI NR), NCBI nucleotide sequences (NT), Kyoto Encyclopedia of Genes and Genomes (KEGG), Swiss-prot, Gene Ontology (GO), euKaryotic Ortholog Groups (KOG) and Protein family (Pfam) were used to de novo annotate all unigenes. The reads of all libraries were deposited in the NCBI SRA (Short Read Archive) with BioProject accession number of PRJNA878928.

Identification and analysis of metabolizing enzyme genes
Identification of genes encoding metabolizing enzymes, CYPs, CCEs and GSTs, of A. custos was firstly performed via keyword searches of the A. custos transcriptome annotation table. The annotated CYP, CCE and GST protein sequences of three bugs, Halyomorpha halys [26], Rhodnius prolixus [27] and Cyrtorhinus lividipennis [15], were used for BLAST queries against the A. custos transcriptome to identify any missed genes. Putative metabolizing enzyme genes were in turn used as queries for additional blasting until no new candidates could be found.
Phylogenetic trees and NCBI blast searches were used to identify subfamilies of A. custos metabolizing enzyme genes. After manually purging small gene fragments, phylogenetic trees were constructed based on the amino acid sequence alignment of metabolizing enzyme genes from A. custos and other bug species by applying MEGA version 6 software [28], using the maximum likelihood method, the Jones-Taylor-Thornton (JTT) model and a bootstrap analysis with 1000 replicates. Subsequently, subfamily assignment for small gene fragments without involving in the phylogenetic trees were performed by the blast searches against NCBI databases and the annotated CYP, CCE and GST protein sequences of H. halys [26], R. prolixus [27] and C. lividipennis [15].

Quantitative real-time PCR (qRT-PCR) validation of RNA-seq data
In order to validate the RNA-seq data, fifteen genes were selected for qRT-PCR analysis using SYBR Premix Ex Taq (Takara Biotechnology Corporation Co. Ltd, Dalian, China) and an ABI Prism 7300 (Applied Biosystems, Foster City, CA). A list of qRT-PCR primers was presented in Table S1. EF1A (Cluster-13084.11009) was selected for the candidate reference genes. Regular PCR assays and melting curve anlyses were performed to verify the specificity of primers ( Figure S1). The qRT-PCR assays were performed with three independent biological replicates and three technical repeats. Relative expression levels for each gene were calculated using the 2 − ΔΔCt method [29]. The correlation analysis of fold change (FC) data (log2 transformation) between RNA-Seq and qRT-PCR was performed using Pearson correlation analysis by SAS (version 8.01, SAS Inc., Cary, USA).

RNA sequencing and transcript assembly
A total of 15 cDNA libraries, three each from guts (G), salivary glands (SG), antennae (A), legs (L) and heads (H), were produced from fifth-instar larvae of A. custos. The raw sequence (with more than 19.4 million reads in each sample) was processed for trimming and filtering, resulting in the Q20 value higher than 97% in each sample (Table S2). The transcriptomics from tissue samples resulted in a total of 41,845 unigenes with a mean length of 1219 bp, and a N50 length of 2183 bp. About 50% of 41,845 unigenes matched known proteins, after annotated with NR, NT, KEGG, Swiss-prot, GO, KOG and Pfam databases. Similarity distribution and species distribution of unigenes that hit in the NCBI NR protein database are shown in Figure S2.

Identification and tissue-specific expression analysis of CYPs
The insect CYP superfamily includes four clans: CYP2, CYP3, CYP4 and the mitochondrial CYP clan (CYPmito), which in turn are subdivided into families and subfamilies [30]. In the A. custos transcriptome, 91 putative genes or gene fragments were annotated as CYPs (Table S3). Based on phylogeny trees and the BLAST searches against NCBI databases and the named CYPs of H. halys [26], R. prolixus [27] and C. lividipennis [15], they were identified and sorted into four CYP clans and 27 families. Specifically, 48 genes were classified to the CYP3 clan, 31 to the CYP4 clan, 5 to the CYP2 clan and 7 to the CYPmito clan ( Table 1). The number of A. custos CYP genes is much less than those from H. halys, R. prolixus and Triatoma infestans, but much more than those from Orius laevigatus, C. lividipennis, Cimex lectularius and Murgantia histrionica (Table 1). A. custos seems to have a larger species-specific expansion on CYP3 and CYP4 clans compared to the other two predatory bugs, O. laevigatus and C. lividipennis ( Table 1).
In order to identify candidated degrading genes responding to xenobiotics and odors, expression patterns among guts, antennae and other tissues/parts of A. custos nymphs were examined based on their FPKM values. To validate the DEGs identified through RNA-seq sequencing, fifteen genes were selected for qRT-PCR analysis. The results of linear regression analysis indicated a significant correlation (R 2 = 0.70, F = 64.35, P < 0.0001) between the data of qRT-PCR and RNA-Seq ( Figure S3).
In Tables 2-4, FPKM values of CYPs in guts, antennae and heads were exhibited. Amongst the 48 CYP3 clan genes, 30 showed significantly higher expression levels in guts than heads, whereas 4 were significantly up-regulated in anntenna than heads, based on a criterion of Log2 (fold) > 1 and P_adj<0.05 for significant difference ( Table 2). Of the 31 CYP4 genes, 7 and 2 were significantly more expressed in guts and anntenna than heads, respectively ( Table 3). All of the 14 members in CYP2 and CYPmito clans exhibited predominant expression levels in heads (Table 4).

Identification and tissue-specific expression analysis of CCEs
Insect CCEs is typically divided into three functional classes: dietary/detoxification (DD), hormone and pheromone processing (HPP) and neuro/developmental (ND) [20,[34][35][36]. A total of 50 putative CCE genes or gene fragments, after annotating the transcriptome and removing a couple of small redundancies (Table S3; Tables 5-6). Our results revealed that DD CCEs are absent in A. custos, which is consistent with recent studies that DD class is absent in hematophagous bugs, predatory O. laevigatus, and herbivory Pentatomidae bugs [31]. Numbers of A. custos CCE genes in the HPP and ND classes are 32 and 18, respectively (Table 5). In the maximum likelihood tree of CCEs, 39 full or nearly full-length sequences of CCEs in A. custos were involved (Fig. 4). The HPP class in A. custos includes secreted β-esterases and juvenile hormone esterases (JHEs). It shows an expansion of the HPP class in A. custos compared with other two predatory bugs ( Table 5). The ND class in A. custos contains catalytic acetylcholinesterases (AchE-1 and AchE-2) and non-catalytic neuroligin, glutactin, gliotactin, neurotactin and uncharacterized clade I proteins (Table 6).
A total of 14 and 7 CCEs were significantly enriched in guts and antennae compared to heads (Table 6). Amongst the 32 HPP class genes, 14 exhibited higher expression levels in guts, whereas 5 were more expressed in antennae than heads (Table 6). Amongst the 15 Fig. 1. Phylogenetic tree of the CYP3 clan genes of A. custos (Ac), Halyomorpha halys (Hh) [26], Rhodnius prolixus (Rp) [27,32] and Cyrtorhinus lividipennis (Cl) [15], constructed by the Maximum likelihood method with the JTT model. ND class genes, almost all of them exhibited higher expression levels in heads than guts or antennae, except for two genes that were more expressed in antennae than heads (Table 6).

Identification and tissue-specific expression analysis of GSTs
Insect GSTs can be grouped into seven classes named Delta, Epsilon, Omega, Sigma, Theta, Zeta, and microsomal [37]. In A. custos, 23 putative GST genes or gene fragments were found (Table S3; Tables 7-8). A phylogenetic analysis of GSTs showed that A. custos GSTs belong to six classes: delta (2 members), omega (2), sigma (12), theta (2), zeta (1) and microsomal (4) (Fig. 5). No epsilon GST gene was found in A. custos. The number of GST genes of different classes found in A. custos mirrors that reported in O. laevigatus (Table 7). Particularly, the sigma class was the largest group in A. custos and other bug species (Table 7).
Compared to heads, 8 out of 23 GSTs showed higher expression levels in guts, including half (6/12) of GST sigma class genes (Table 8). Only two GSTs were significantly more expressed in antennae than heads (Table 8).

Discussion
In the present study, 91 CYPs, 50 CCEs and 23 GSTs were identified in a transcriptome analysis of different tissues or parts from fifth-instar nymph of A. custos. Notably, almost all of these genes of A. custos were more identical to their orthologues of H. halys in comparsion with those of C. lividipennis, R. prolixus and other bugs, supporting that predatory stink bugs were evolved from the phytophagous stink bugs [38,39]. The number of detoxification famliy gene members varying across insect species might be associated with the functional differentiation during insect evolution and related to the process of adaptation to the environment [40]. After comparing the number of these three famliy gene sequences in A. custos with those of the corresponding genes in other several Heteroptera species, we found a similar pattern of the member numbers of CYPs, CYP3, CYP4 and the HPP CCE class across Heteroptera species. A. custos contained more detoxification clan or class genes than other predatory bugs, O. laevigatus and C. lividipennis, but fewer than the herbivory stink bug H. halys and the hematophagous bug R. prolixus. The herbivory stink bug H. halys embracing the more expanded and diverse set of detoxification genes might due to its extreme generalist behavior [26,41]. Although A. custos, O. laevigatus and C. lividipennis are zoophytophagous that are carnivores and sometimes feed on plants, A. custos might need more detoxification genes to adapt more habitats and feed on more prey species [42]. These results indicate that CYP3, CYP4 and HPP CCEs relatively Note: * indicated that the expression level of one gene in guts (G) or antennae (A) was significantly higher than that in heads (H), with log2(fold) > 1 and P_adj<0.05. The heatmap scale was based on log10 (FPKM+1) values:Min = 0.0 Max = 2.47.

Fig. 2.
Phylogenetic tree of the CYP4 clan genes of A. custos (Ac), Halyomorpha halys (Hh) [26], Rhodnius prolixus (Rp) [27,32] and Cyrtorhinus lividipennis (Cl) [15], constructed by the Maximum likelihood method with the JTT model.  expanded in A. custos might be importance detoxification gene groups responding to more complex environment. Our results also showed the different tissue-specific expression patterns of CYPs, CCEs and GSTs in A. custos, which provide useful information for the identification of candidated degrading genes responding to xenobiotics and odors. The insect gut serves as a major digestion organ, involved in digestion and detoxification of food and xenobiotics. A total of 37 CYPs, 14 CCEs and 8 GSTs were enriched in guts of A. custos, suggesting that they were proposed to have important roles on xenobiotic detoxification. In addition, 6 CYPs, 5 CCEs and 2 GSTs were up-regulated in A. custos antennae, indicating that these antenna-enriched metabolizing genes might have the role on ordorant degradation.
CYPs can readily catalyze the detoxification of chemicals by hydroxylation, epoxidation, dealkylations and a great variety of other monooxygenase reactions [43]. Our finding that most of CYP3 genes (37/48) were mainly expressed in guts of A. custos nymphs, showing that CYP3 clan members were the major detoxification gene group in the predatory stink bugs. Expression levels of CYP4 clan members were quite diverse among different tissues of A. custos nymphs. Seven gut-specific CYP4 genes of A. custos might be related to xenobiotic detoxification, whereas two CYP4 genes enriched in A. custos antennae might be ODEs.
The well-studied type of ODEs are P450s. The antennal enriched CYPs not only degrade plant volatiles and insecticides but also inactivate the pheromones in olfactory organs [13]. An antennal specific enzyme Dendroctonus ponderosae CYP345E2 that was the first olfactory CYP functionally characterized in vitro in an insect may catalyze the oxidation of many volatile monoterpenes for odorant clearance [44]. An antennal CYP gene Spodoptera litura CYP4L4 played a role in the recognition sex pheromones [45]. The   Fig. 3. Phylogenetic tree of the CYP2 clan and CYPmito genes of A. custos (Ac), Halyomorpha halys (Hh) [26], Rhodnius prolixus (Rp) [27,32] and Cyrtorhinus lividipennis (Cl) [15], constructed by the maximum likelihood method with the JTT model.

Table 5
Numbers of CCE genes annotated in A. custos (in this study), Orius laevigatus [31], Cyrtorhinus lividipennis [15], Rhodnius prolixus [27,32], and Halyomorpha halys [26]. *The number of CCE genes in Halyomorpha halys was revised to 76 after deleting 6 alternative splicing variants based on Sparks et al. [26]. identification of ten antennal specific CYPs of A. custos in the present study will be informative in the future studies about the chemical ecology of the predatory stink bugs. Generally, insect CYPs in the CYP2 and mito clans have been regarded as conserved players in cell signaling and developmental processes, which is consistent with our observation that they were mainly expressed in heads of A. custos nymphs. The CYP2 and CYPmito clans contain several halloween genes for ecdysone and 20-hydroxyecdysone (20-E) synthesis (CYP302A1, CYP306A1, CYP307B1, CYP314A1 and CYP315A1), CYP18 for 20-E turnover, and CYP15 for juvenile hormone synthesis, CYP301 for cuticle formation, CYP303 gene for the structure and function of sensory organs, and CYP305 for unknown function [26]. But a recent study confirmed that some of them are also capable of metabolizing xenobiotics [46]. The potential roles of CYP2 and mito clan P450s of predaroty bugs in insect chemical defense remain unclear.
Recent studies suggested that DD CCE class is absent in hematophagous bugs, predatory O. laevigatus, and herbivory Pentatomidae bugs [31]. Our results revealed that DD CCEs are also absent in the predatory Pentatomidae bug A. custos, suggesting the lose of DD  CCEs during the evolution of Heteroptera. In A. custos, most of known ND CCE class genes exhibited higher expression levels in heads, which suggests a potential function in neural development and signal transduction for these genes. Almost half of HPP CCE clade members (14/32) showed higher expression levels in guts of A. custos, suggesting that these HPP CCE genes might be involved in detoxification. What's more, five HPP CCE up-regulated in antennae might be ODEs. CCEs is also one of the most well documented ODEs [47]. The first CCE enzyme related to olfaction in insects was an antennae-specific esterase (ApolPDE) from Giant silk moth, Antheraea polyphemus [48]. Several more antennal CCEs involved in odorant degradation in insects have been identified subsequently, such as esterase 6 (EST-6) in Drosophila melanogaster [49,50], SlitCXE7 and SlitCXE10 in Spodoptera littoralis [51,52], SexiCXE4 and SexiCXE14 in S. exigua [53,54]. The expression in a specific tissue could shed lights on the function of a CCE. Several antennal-specific CCEs in A. custos might be involved in the hydrolysis of host  [26] and Rhodnius prolixus (Rp) [27,32] constructed by the Maximum likelihood method with the JTT model. Two Drosophila melanogaster (Dm) dietary/detoxification (DD) class genes and four Nilaparvata lugens (Nl) DD class genes were added into the tree, in order to showing the absent of DD class in Heteroptera. Expasion of the hormone and pheromone processing (HPP) class in bugs were shown in the tree.
GSTs that are crucial in odorant degradation as the olfactory genes of insects generally show preferential expression in the antennae, for example, GST-msolf1 of M. sexta [64], GST-pxcs1 of P. xuthus [65], GmolGSTD1 of Grapholita molesta [66] and SzeaGSTd1 of Sitophilus zeamais [67]. In A. custos, only two antennae-specific GSTs were identified, whose functions related to odorant degradation needed be further investigated.
Predatory stink bugs are zoophytophagous generalists of various insect pests in agriculture and forest plantations [68]. A. custos feeds on many species of pest insects and perfers to live in elm, poplar, cotton, soybean and other plants [2]. Stink bugs (Pentatomidae) produce unpleasant scent odors that function as alarm and defense signals against natural enemies [69]. Therefore the predatory stink bugs depend on their antennae to perceive a diversity of airborne chemical cues, including prey odorants, plant volatiles, sex pheromones and scents from stink bugs, to find palatable preys, preferable host plants and conspecific partner, and to avoid interspecific competition and natural enemies. These insects rapidly respond to stimuli in their environment, mainly relying on the termination of odorant signals from the olfactory sensilla following stimulation by ODEs [13]. A number of candidate A. custos ODEs that were identified in the present study will be informative in the future studies about the elucidation of their olfactory system and the improvement of their bio-control efficiency.

Conclusions
In the present study, the comprehensive analysis of xenobiotics detoxification and odor degradation genes in the predatory stink bug were conducted. A total of 91 CYPs, 50 CCEs and 23 GSTs were identified based on RNA-seq transcriptomic analysis. Tissuespecific expression analysis showed that most of gut-enriched detoxification enzyme genes belonged to CYP3 and CYP4 clans, HPP CCE and GST Sigma classes, whereas antennae-specific odorant degrading genes were contained in CYP4 clan and CCE HPP class. A further in-depth study examined these xenobiotics detoxification and odor degradation genes will provide basic information for studying the improvement of the bio-control efficiency of the predatory stink bugs.

Author contribution statement
Wenhong Li: Conceived and designed the experiments; Performed the experiments; Wrote the paper.

Data availability statement
Data associated with this study has been deposited at NCBI SRA with BioProject accession number of PRJNA878928.

Additional information
Supplementary content related to this article has been published online at [URL].

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  [26], Rhodnius prolixus (Rp) [27,32] and Cyrtorhinus lividipennis (Cl) [15], constructed by the Maximum likelihood method with the JTT model. Expasion of the Sigma class in bugs were shown in the tree.